#include "atmoslayer.h"
#include <math.h>
#include <QTime>
#include <QDebug>

//#define INTENSE_DEBUG

/***************************************************
 for lightUp: lightUp[0][Azi] is zenith 0
 for lightD:  lightD[0][Azi] is zenith Pi(Pi-0),
	      lightD[Pi/2][Azi] is zenith Pi/2
 ***************************************************/


void AtmosLayer::solarScatCoeff(int iThetaS, int iPhiS,
				linearPolar *lightOut,
				int scatMode)
{
}

//value tested against Mathematica, 2-8
void AtmosLayer::scatCoeff(int iThetaS, int iPhiS,
			   int iThetaI, int iPhiI,
			   double lightIncident[3],
			   linearPolar *lightOut, int scatMode)
{
}

//#define SCAT_CALL_DEBUG
void AtmosLayer::upInUpScat()
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "AtmosLayer::upInUpScat()";
#endif
	int iThetaS, iPhiS, iThetaI, iPhiI;
	double muI, muS;
	for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){

	    muS = Cos(iThetaS*M_PI/N_PI + ANGLE_STEP/2);

	    for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

		m_lightScatUp[iThetaS][iPhiS][I]+=
			m_lightInUp[iThetaS][iPhiS][I] *
			exp(-m_tau/muS);
		m_lightScatUp[iThetaS][iPhiS][Q]+=
			m_lightInUp[iThetaS][iPhiS][Q] *
			exp(-m_tau/muS);
		m_lightScatUp[iThetaS][iPhiS][U]+=
			m_lightInUp[iThetaS][iPhiS][U] *
			exp(-m_tau/muS);

		for(iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){

		    muI = Cos(iThetaI*M_PI/N_PI + ANGLE_STEP/2);
		    double coeffPath;
		    coeffPath = muI/muS *
			    (1 - exp(-1 * m_tau / muI))*
			    Sin(iThetaI*M_PI/N_PI + ANGLE_STEP/2)*
			    pow(M_PI/N_PI,2);

		    for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){

			linearPolar coeff;

			scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInUp[iThetaI][iPhiI], &coeff, UP_IN_UP_SCAT);

			m_lightScatUp[iThetaS][iPhiS][I] +=
				coeffPath * coeff.i ;

			m_lightScatUp[iThetaS][iPhiS][Q] +=
				coeffPath * coeff.q;

			m_lightScatUp[iThetaS][iPhiS][U] +=
				coeffPath * coeff.u;

		    }
		}
	    }
	}
}

void AtmosLayer::upInDownScat()
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "AtmosLayer::upInDownScat()";
#endif
    int iThetaS, iPhiS, iThetaI, iPhiI;
    double muI, muS;
    for(iThetaS = 0; iThetaS > -N_ZENITH; iThetaS--){

	muS = qAbs(Cos(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2)); //Pi - theta

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    for(iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){

		muI = Cos(iThetaI*M_PI/N_PI + ANGLE_STEP/2);
		double coeffPath;
		coeffPath = muI/muS *
			(1 - exp(-1 * m_tau / muI))*
			Sin(iThetaI*M_PI/N_PI + ANGLE_STEP/2)*
			pow(M_PI/N_PI,2);

		for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){

		    linearPolar coeff;

		    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInUp[iThetaI][iPhiI], &coeff, UP_IN_D_SCAT);

		    m_lightScatD[-iThetaS][iPhiS][I] +=
			    coeffPath * coeff.i ;

		    m_lightScatD[-iThetaS][iPhiS][Q] +=
			    coeffPath * coeff.q;

		    m_lightScatD[-iThetaS][iPhiS][U] +=
			    coeffPath * coeff.u;
		}
	    }
	}
    }
}

//downward scatter-flux response caculate
void AtmosLayer::downInUpScat()
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "AtmosLayer::downInUpScat()";
#endif
    int iThetaS, iPhiS, iThetaI, iPhiI;
    double muI, muS;
    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	muS = Cos(iThetaS*M_PI/N_PI + ANGLE_STEP/2);

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    for(iThetaI = 0; iThetaI > -N_ZENITH; iThetaI--){

		muI = qAbs(Cos(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2));
		double coeffPath;
		coeffPath = muI/muS *
			(1 - exp(-1 * m_tau / muI))*
			Sin(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2)*
			pow(M_PI/N_PI,2);

		for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){

		    linearPolar coeff;

		    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInD[-iThetaI][iPhiI], &coeff, D_IN_UP_SCAT);

//		    qDebug() << "before" <<  iThetaI << iPhiI << iThetaS << iPhiS << m_lightScatUp[iThetaS][iPhiS][I] << coeff.i << coeff.q << coeff.u;

		    m_lightScatUp[iThetaS][iPhiS][I] +=
			    coeffPath * coeff.i ;

		    m_lightScatUp[iThetaS][iPhiS][Q] +=
			    coeffPath * coeff.q;

		    m_lightScatUp[iThetaS][iPhiS][U] +=
			    coeffPath * coeff.u;

//		    qDebug() << "end   " <<  iThetaI << iPhiI << iThetaS << iPhiS << m_lightScatUp[iThetaS][iPhiS][I];

		}
	    }
	}
    }
}

//downward scatter-flux response caculate,
void AtmosLayer::downInDownScat()
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "AtmosLayer::downInDownScat()";
#endif
    int iThetaS, iPhiS, iThetaI, iPhiI;
    double muI, muS;
    for(iThetaS = 0; iThetaS > -N_ZENITH; iThetaS--){

	muS = qAbs(Cos(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2));

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    m_lightScatD[-iThetaS][iPhiS][I]+=
		    m_lightInD[-iThetaS][iPhiS][I] *
		    exp(-m_tau/muS);
	    m_lightScatD[-iThetaS][iPhiS][Q]+=
		    m_lightInD[-iThetaS][iPhiS][Q] *
		    exp(-m_tau/muS);
	    m_lightScatD[-iThetaS][iPhiS][U]+=
		    m_lightInD[-iThetaS][iPhiS][U] *
		    exp(-m_tau/muS);

	    for(iThetaI = 0; iThetaI > -N_ZENITH; iThetaI--){

		muI = qAbs(Cos(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2));
		double coeffPath;
		coeffPath = muI/muS *
			(1 - exp(-1 * m_tau / muI))*
			Sin(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2)*
			pow(M_PI/N_PI,2);

		for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){

		    linearPolar coeff;

		    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInD[-iThetaI][iPhiI], &coeff, D_IN_D_SCAT);

		    m_lightScatD[-iThetaS][iPhiS][I] +=
			    coeffPath * coeff.i ;

		    m_lightScatD[-iThetaS][iPhiS][Q] +=
			    coeffPath * coeff.q;

		    m_lightScatD[-iThetaS][iPhiS][U] +=
			    coeffPath * coeff.u;

//		    if(iThetaS == -10 && iPhiS == 10)
//			qDebug() << iThetaI << iThetaS << iPhiI << iPhiS << m_lightScatD[-iThetaS][iPhiS][I] << coeff.i;
		}
	    }
	}
    }
}



//#define SOLAR_SCAT_DEBUG
void AtmosLayer::solarUpScatCaculate()
{
    if(m_solarIntensity == 0)
	return;

    int iThetaS, iPhiS;
    double muI, muS;

    muI = qAbs(Cos(m_solarZen));

    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){

	muS = Cos(iThetaS*M_PI/N_PI + ANGLE_STEP/2);

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    linearPolar coeff;

	    solarScatCoeff(iThetaS, iPhiS, &coeff, D_IN_UP_SCAT);

	    m_lightScatUp[iThetaS][iPhiS][I] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI)) *
		    coeff.i ;

	    m_lightScatUp[iThetaS][iPhiS][Q] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    coeff.q ;

	    m_lightScatUp[iThetaS][iPhiS][U] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    coeff.u ;


	}
    }

}

void AtmosLayer::solarDownScatCaculate()
{
    if(m_solarIntensity == 0)
	return;

    int iThetaS, iPhiS;
    double muI, muS;

    muI = qAbs(Cos(m_solarZen));

    for(iThetaS = 0; iThetaS > -N_ZENITH; iThetaS--){

	muS = qAbs(Cos(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2));

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    linearPolar coeff;

	    solarScatCoeff(iThetaS, iPhiS, &coeff, D_IN_D_SCAT);

	    m_lightScatD[-iThetaS][iPhiS][I] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    coeff.i ;

	    m_lightScatD[-iThetaS][iPhiS][Q] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    coeff.q ;

	    m_lightScatD[-iThetaS][iPhiS][U] +=
		    muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    coeff.u ;

	}
    }

}
